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ABSTRACT. We review recent achievements in the solution of the initial- 
value problem for quantum back- reaction in scalar and spinor QED. The 
problem is formulated and solved in the semiclassical mean-field approxima- 
tion for a homogeneous, time-dependent electric field. Our primary motiva- 
tion in examining back-reaction has to do with applications to theoretical 
models of production of the quark-gluon plasma, though we here address 
practicable solutions for back-reaction in general. We review the application 
of the method of adiabatic regularization to the Klein-Gordon and Dirac 
fields in order to renormalize the expectation value of the current and derive 
a finite coupled set of ordinary differential equations for the time evolution of 
the system. Three time scales are involved in the problem and therefore cau- 
tion is needed to achieve numerical stability for this system. Several physical 
features, like plasma oscillations and plateaus in the current, appear in the so- 
lution. From the plateau of the electric current one can estimate the number 
of pairs before the onset of plasma oscillations, while the plasma oscillations 
themselves yield the number of particles from the plasma frequency. 

We compare the field-theory solution to a simple model based on a rela- 
tivistic Boltzmann-Vlasov equation, with a particle production source term 
inferred from the Schwinger particle creation rate and a Pauli-blocking (or 
Bose-enhancement) factor. This model reproduces very well the time behav- 
ior of the electric field and the creation rate of charged pairs of the semiclassi- 
cal calculation. It therefore provides a simple intuitive understanding of the 
nature of the solution since nearly all the physical features can be expressed 
in terms of the classical distribution function. 
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1 Introduction 



This review deals with the physical situation in which initially there exists 
a classical electric field out of which pairs of charged particles tunnel. The 
particles are accelerated by the field, producing a current and, in the sequel, 
a field that acts counter to the original one. Eventually, plasma oscillations 
are set up. This so-called back-reaction problem has been of considerable 
interest in recent years in particle physics, serving for example as a model 
of the production of the quark-gluon plasma. One supposes that when two 
nuclei collide at ultra-relativistic energies, they induce color charges on each 
other. Passing through each other, they leave in their wake a color electric 
field from which quarks and gluons emerge more or less according to the above 
scenario. Of course, the color electric field is non-Abelian, whereas we take it 
to be Abelian in the version of the back-reaction problem we defined above. 
Furthermore, the field does not fill all space homogeneously, a condition that 
thus far one has had to impose in order to carry out practical calculations, 
as we shall discuss below. 

Back reaction has been studied intensively in recent years in the domain of 
inflationary cosmology, where instead of a time-varying electric field one has 
the time-dependent gravitational metric. Again pairs emerge via tunneling 
and then act back on the field, in this case through their masses. 

Our own interest in the back-reaction problem was sparked by its applica- 
tion to the quark-gluon plasma problem, and it is therefore almost exclusively 
in these terms that we shall present our discussion here. Much of the evo- 
lution of regularization and renormalization methods for the back-reaction 
problem [1-17] has taken place, however, in the context of inflationary cos- 
mology. In this area there exists a sizable literature, summarized — except, 
of course, for the more recent developments — in the books of Birrell and 
Davies [18] and of Fulling [19]. Earlier surveys of the status of this problem 
can be found in the reviews of Zel'dovich [20] and of DeWitt [21]. We shall 
not describe these developments here, but rather pick up our discussion with 
the treatment of adiabatic regularization given by Cooper and Mottola [16], 
which has proved an adequate starting point for the practicable calculational 
schemes upon which we focus. 

In concentrating on methods that lead to workable calculations for back- 
reaction, we shall not address the vast number of papers that deal primarily 
or exclusively with the Schwinger mechanism [22-27] and its many applica- 
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tions such as the flux-tube model 2 [29-32]; nor shall we touch on the gen- 
eralizations of this mechanism to handle pair production in given, external, 
time-dependent electric fields, reviewed in Refs. [33, 34]; nor do we deal with 
the various efforts to consider pair production in restricted volumes [35]. We 
shall restrict ourselves to methods based on adiabatic regularization for the 
renormalization of the electric current or on methods that have been shown 
[5, 9, 18, 19] to be equivalent to it. 3 

Interest in back-reaction in the context of the quark-gluon plasma 4 (QGP) 
arose through work of Bialas and Czyz and coworkers [39, 40], of Kajantie 
and Matsui [41], and of Gatoff, Kerman, and Matsui [42]. In these series 
of studies, the authors modeled the production and early dynamics of the 
QGP by considering a transport equation for quarks or gluons, on the right- 
hand side of which appears a source term based on the Schwinger formula 
for pair production in an electric field fixed in time. The electric field is, 
in turn, governed by Maxwell's equations, with the current of the produced 
pairs figuring in the time evolution of the electric field. This coupled system 
is solved, and the subsequent development exhibits plasma oscillations and 
other distinctive behavior (as we shall see in our discussion here). 

The genesis of the classical transport equation, while intuitively very ap- 
pealing, is far from self-evident: The whole point in back-reaction is that 
the electric field must change in response to the current of produced pairs, 
so that the assumption of a fixed field, inherent in the Schwinger formula, 
is out of place. Furthermore, when one imagines deriving a transport equa- 
tion directly from the field equations through the well-known procedure of 
the Wigner transformation [43], it is clear that the equation cannot lead di- 
rectly to a source term. 5 Since the Wigner treatment is completely general, 
it follows that in any quantum mechanical treatment there is no possible 
separation between the role of the electric field in accelerating charged par- 
ticles and its role in producing these particles out of the vacuum through 
tunneling. Thus it is quite mysterious how a classical transport equation can 
arise in this formalism. 

In order to examine back-reaction in detail, we return to the original 

2 See also [28]. 

3 See also the early calculation in one spatial dimension by Ambj0rn and Wolfram [36]. 
4 For general introductions to this subject, see, e.g., [37, 38] 

5 See also Refs. [44-49] for discussion of this point, which lies beyond the scope of the 
present review. 
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initial-value problem in quantum field theory. We write a coupled set of 
field equations for charged particles in interaction with an electric field. For 
purposes of comparing with the models [39-42] for the QGP, it is adequate 
to take this electric field as classical. (Purists may note that this is the 
leading order in a 1/N expansion, where N is the number of flavors of the 
matter field.) Clearly one will wish eventually to consider a quantized elec- 
tromagnetic field (i.e., higher orders in 1/N) so as to incorporate the essential 
physics of radiation for more realistic applications. 6 The specification of ini- 
tial conditions for the matter field reduces the field equations to c-number 
differential equations. This development was carried out initially for bosons 
[16, 50], and is presented here in Section 2. 

The numerical solution of these coupled equations is confronted by two 
obstacles: (i) The field equations must be renormalized in order to render 
them finite, and (ii) it emerges that the calculated momentum distribution 
of the produced particles is highly oscillatory so that great care is required 
to produce reliable numerical results. The first of these problems is solved by 
taking as a starting point the extensive work noted earlier in renormalizing 
the pair-production problem for inflationary cosmology [1-19], though, as 
we shall see, some modifications to this procedure are required in practice. 
Unfortunately, the available techniques limit us to the case of a spatially 
homogeneous electric field. Overcoming the second obstacle merely requires 
sufficient computer power to achieve long-range numerical stability in spite of 
short-range oscillations. These points are developed in some detail in Section 
2. We there review results [50] in 1 + 1 dimensions and present results in 
3 + 1 dimensions for the first time. 

Somewhat to our surprise, the classical transport equation used [39-42] 
for the QGP permits very nearly a quantitative simulation of the field-theory 
result, thus considerably bolstering one's belief in the legitimacy of various 
detailed steps involved at the technical level in both procedures. The agree- 
ment between the two procedures is greatly improved by including in the 
transport equation a term that provides for Bose-Einstein enhancement in 
the production process; this term is found to have a particular form which is 
motivated by careful consideration of pair production in an external classical 
electric field. 

6 There is also the still more difficult issue of the non-Abelian field involved in the QGP; 
this is also of higher order in 1/N. 
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The treatment of fermion pair production [51], presented in Section 3, 
involves some technical issues concerning the use of the adiabatic regulariza- 
tion procedure. The physics of Pauli blocking makes it rather less evident 
that a classical transport equation can succeed here, so it is, perhaps, even 
more surprising than for bosons that a classical transport equation gives an 
adequate description of the results, even at a quantitative level. In this case, 
good agreement requires the introduction of an explicit Pauli-blocking factor 
in the transport equation. 

The application to the QGP requires going beyond the case of a uniform 
electric field. Luckily, one can take advantage of the approximate invariance 
under longitudinal boosts which is often assumed in modeling particle pro- 
duction in the central rapidity region [52, 53]. The transformation of the 
field equations to comoving coordinates allows renormalization just as in the 
case of the homogeneous electric field, and numerical solution follows [54]. 
We do not review these developments here. 

The work covered here treats a highly-idealized situation of a classical 
electric field, homogeneous in space and producing particles that possess no 
mutual interactions. This can at best serve as a starting point for more 
serious studies that contain the complications that are essential for physical 
applications. However, we feel that the beginning phase is now more or less 
complete, and therefore is deserving of review here. The incorporation of 
particle scattering and of bremsstrahlung is straightforward, and such work 
is in progress. Lifting the restriction of spatial homogeneity (or of boost 
invariance), on the other hand, appears very difficult. 

2 Boson pair production 
2.1 Introduction 

Beginning with Sauter's pioneering paper [22] in 1931, the problem of spon- 
taneous pair production in the presence of an external electric field has been 
investigated by many authors [23-27, 33, 34]. The most commonly used 
formula for the spontaneous pair creation rate per unit volume, derived by 
Schwinger [24] , is based on a field-theory calculation for a constant and ho- 
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mogeneous electric field E, to all orders in E, 

-P' +1 Wy E|exp(-^ / B ), (2.1) 

where s is the spin of the particles produced, (3 n = (— l)™ -1 for bosons and 
(3 n — 1 for fermions, and E = m 2 /e. This formula has been used extensively 
in modeling particle production in the central rapidity region in high-energy 
nucleus-nucleus collisions [39-42, 55, 56]. 

In this expression the electric field is held constant by an external agent. 
Thus, there is no effect of the produced particles on the original electric 
field — no back-reaction. Moreover, the mutual interactions of the particles 
are not taken into account. As long as the intensity of the initial electric field 
is small, E<^Eq, Schwinger's formula is valid. However, when it is applied for 
a strong electric field, E > Eq, it gives w oc (E/E0) 2 , violating unitarity for 
sufficiently strong fields [33]. In such circumstances it is necessary to include 
both the interaction between the particles created and the screening of the 
external field. 

As discussed in Section 1, we assume that the electric field is classical, 
Abelian, and homogeneous. We impose the restriction to a spatially uniform 
electric field not just for simplicity, but also to dispose of infinities that would 
otherwise appear in the expectation value of the current operator in the 
Maxwell equation. These infinities are removed by adiabatic regularization 
and renormalization. It is important to note that in an initial value problem 
the divergences in the expectation values of the electric current or the energy- 
momentum tensor are time dependent. Nonetheless, the charge and mass 
renormalizations are time-independent and the infinities which appear in 
the currents are products of time-independent infinite quantities with time- 
dependent finite quantities. 

A scheme for solving the quantum back-reaction problem in scalar QED 
was put forth by Cooper and Mottola [16]. Some analytical progress based 
on this scheme has been made [17], but is limited to very short times or 
very weak electric fields. In order to see substantial pair production and the 
subsequent onset of plasma oscillations we need to begin with a strong initial 
electric field, E ~ Eq, and to follow the evolution to large times, t^m/eE. 
It is useful to investigate the problem in 1+1 dimensions as a first step in 
testing numerical procedures since there are no transverse momenta to take 
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into account, and since renormalization is relatively trivial. We then proceed 
to 3 + 1 dimensions. 

Adiabatic regularization will only render the theory finite if we impose a 
special initial configuration of the charged matter field — the adiabatic vac- 
uum. This does not mean that there are no particles in the initial state. As 
long as the matter field and the electric field interact with each other, the 
expectation value of the number operator is not equal to the asymptotic par- 
ticle number density, and even this adiabatic vacuum may contain a nonzero 
density. Nonadiabatic initial conditions would introduce infinities into the 
time derivatives of the current at t — 0. Ultimately, the physical behavior 
we find in our calculations offers us a way to determine an effective particle 
density during the evolution of the system 7 . 

2.2 Pair production in semiclassical scalar QED 
Equations of motion and second quantization 

The scalar QED coupled equations of motion in the mean-field approximation 
(and without a 4 self- interaction term as included in [16]) are given by the 
Klein-Gordon equation 

[(9" + ieA»){d» + ieAfj) + m 2 ]$(x) = (2.2) 

and the semiclassical Maxwell equations 

d,F^ = (i\f\i), (2.3) 

where \i) is the initial configuration of the charged matter field, and the 
current j v of the charged scalar field is symmetrized in order to be odd 
under charge conjugation, 

f = -[$t _ (£>"$)t$ _ $(zr$) f + (d"$)&}. 

Here D u = d u + ieA u is the covariant derivative. Note that the expectation 
value (i\j u \i) does not vanish because \i) cannot be an eigenstate of the 
charge-conjugation operator as long as a time- varying electric field is present. 

7 One may also add explicitly a finite density of particles to the initial adiabatic vacuum 
without disturbing the finiteness of the theory. 
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As we shall see this expectation value is divergent, and available practical 
schemes which remove this divergence are, at the moment, limited to the case 
of spatially homogeneous electric fields [18, 19]. This homogeneity reduces 
the Maxwell equations to a single equation 

= mo- (2.4) 

The scalar charge density j° vanishes everywhere because of Gauss' Law, and 
the magnetic field is a constant which we choose to be zero. We have selected 
the gauge A = 0, and we choose E to lie along the z axis. Thus A z is the 
only nonvanishing component, A u = (0, 0, 0, A(t)). 

Spatial homogeneity enables us to express the boson field operator $ as a 
Fourier integral in which the coefficients are the Heisenberg operators a^(t) 
and b-]J(t), 

t) = J :J= e lk x [a k (f) + 6_ k t(*)] , (2.5) 

where u;£ = \fk? + rn 2 and [dk] = d d k/(2ir) d , with d the spatial dimension. 
Using the Klein-Gordon equation (2.2), the operator O k (t) = a k (i) + &_ k (£) 
satisfies 

^ + ^(t)O k = 0, (2.6) 

where 

u k (t) 2 = [k-eA(t)] 2 + m 2 . (2.7) 

Since u; k (£) is a c- number, the annihilation and creation operators a k (i) and 
b-yj(t) can be expressed in terms of ak(0) and 6_k (0) by means of a Bo- 
golyubov transformation, 

a k (t) = w k (t)a k (0)+t; k (t)6_ k t (0) 
b-J(t) = v k *(t)a k (0)+u k *(t)b^(0). (2.8) 

It is convenient to rewrite <£> in terms of mode amplitudes, 
namely, 

t) = J [dk] [/ k (f)a k (0) + &(t)bU(0)] e lk ' x . (2.9) 
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Thus /k satisfies Eq. (2.6), 

d 2 m 



+ «£(t)/k(t)=0 . (2.10) 



The Bogolyubov coefficients are related to the /k(£) by 

«k(*) = 7== and v k (t) = = — . (2.11) 

The canonical commutation relations 

[$(t,x),n(t,y)]=^(x-y) (2.12) 

imply via (2.5) that 

Mf),at(f)] = [6 k (t),6i(t)] = (27r)^(k- q) . 
These relations will be satisfied if 

M*)| 2 -Mt)| 2 = l , (2.13) 

or, by (2.11), 

/k/k - /k/k = i. (2.14) 

Since the Wronskian of Eq. (2.10) is time independent, this condition is 
satisfied automatically, provided that it is satisfied by the initial conditions 
/k(0),/ k (0). 

Eq. (2.14) implies that fk(t) can be expressed in terms of the real and 
positive function flk(t), 



/k(f) = r=== ex P 

V / 2^ k (t) 



n k (t')dt' 



(2.15) 



if Jlk(0 satisfies 



\2 ^k 3 / f2 k \ 2/ 



The representation (2.15) of / k in a WKB-like form enables us to obtain the 
WKB-like equation (2.16) for fi k , and will help us identify divergences in the 
expectation values of the current and the energy-momentum tensor. 
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The adiabatic expansion and the adiabatic vacuum 



The effect of pair production is related [33] to the appearance of nega- 
tive frequencies in the solution for fk(t). Looked at naively, however, the 
parametrization (2.15) might seem to have only positive frequencies when 
the interaction with the electromagnetic field is switched off. In fact, this is 
not the case since Qk(t) for large t does not approach a constant. To see this, 
we express fk(t) in a different manner using the addition formula [57] 



i 



exp 



cik exp 



+ c 2k exp 



i J 1 n k (t')dt 

-i f W k (t')dt' 



where the constants Cik and c 2 k fulfill 



klkl 2 - |c 2 k| 2 = 1 



(2.17) 



(2.18) 



and where W^(t) satisfies the same equation of motion (2.16) as does fi k (t). 
According to a theorem in asymptotic expansions, which is related to the 
generalized WKB or Liouville-Green approximation, if uj e C°°, the solutions 
of (2.10) are linear combinations of basis solutions [the right-hand side of 
(2.17)] which can be approximated as [19, p. 154] 



1 



exp 



exp 



-i J W k (t')dt' 
Ti fw?{t')dt' 



2W^{t) 

^ k (t) 2 [l + 5 2k (t)Lu k 2 + 5 4k (t)cu k 4 + .. 



(2.19) 



Here <5„k is a function of c^k and its derivatives at t up through ou k n \t), and 
is determined by the order at which the series is truncated. The coefficients 
5 n k are bounded as uo k — > oo. Furthermore, the error term is uniform in t on 
finite intervals. In a static region, when the interaction is switched off (as is 
expected at large times), <5„k = 0. Therefore, it is obvious that the general 
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solution fk(t) of (2.10), where both ci k , C2k 7^ 0, exhibits negative frequencies 
as well as positive ones. 

The amplitude f k (t) can also be represented by a combination of the 
WKB solution and its complex conjugate with time-dependent coefficients 



2W£ d (t) 



a k (t)e-*f w ^ dt ' +p k (t)e +l f 



W£ d (t')dt' 



(2.20) 



where W£ d is the generalized WKB function, up to the desired order. For 
large momentum k, or at times when the electromagnetic interaction is ab- 
sent, we have W k (t) rs W{f(i). In this adiabatic limit W£ d (t) is nearly a 
constant, and a k and f3 k must be constant up to the adiabatic order N to 

which we expanded W£ d (t), because 2W£ d (t) 1/2 exp (=p f W^(t')dt'^ are 
solutions of (2.10) to this order. If we choose a k (t ) = 1 an d Ac(^o) = for 
an initial time t it follows that 

a k (t) = l + 0(iu^ N " 1 ) and (3 k (t) = + Ofa"' 1 ) (2.21) 

for all times. The exact modes f k which correspond to this choice are denoted 
"adiabatic positive-frequency modes of adiabatic order N." Once we choose 
«k(^o) = 1 and /3 k (t ) = 0, then, for large momentum, Cik and C2k are 

Clk = 1 + O^' 1 ) and c 2k = + O^"" 1 ). (2.22) 

Note that at t — to the exact modes are equal to their adiabatic approxima- 
tion of order N. As is seen from the above theorem, the WKB approximation 
for W k is asymptotically correct. 

We see above that the exact adiabatic positive-frequency modes can 
be written as a linear combination (2.17), and for our particular choice, 
«k(^o) = 1 and P k (to) = 0, condition (2.22) is satisfied. An asymptotic ex- 
pansion of modes Q k may therefore differ from the expansion of W k only in 
order O^^^ 1 ). Thus it is legitimate to extract the large k behavior of f2 k 
by adiabatic expansion, even though Q k contains elements of both positive 
and negative frequencies. We also equate the time derivatives of the exact 
modes with the time derivatives of the iVth order approximate adiabatic 
modes at t — to- The annihilation operators a k and b k that correspond to 
these modes (which are matched with their Nth adiabatic approximations 
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at t — to) annihilate a vacuum state, which is said to be the "adiabatic vac- 
uum" \0)a- If is asymptotically static, then f^ 1 and /£ ut are m t ne class 
of solutions approximated by the positive-frequency modes f£ for any N. 
In this special case the adiabatic vacuum can be identified with the usual 
vacua |0 in ) and |0 out ), where the exact modes are matched to their adiabatic 
approximation at t — > — oo and t — > +00, respectively. 

As we shall see, our initial conditions for the back-reaction problem are 
not static because of the presence of an electric field at t — t . It is important 
to note that this adiabatic vacuum is not precisely a vacuum state, that is, 
an inertial particle detector would register a bath of quanta when the field is 
in this adiabatic vacuum. Furthermore, the matching of the mode functions 
may be specified at any time to, and therefore there is no unique iV th order 
adiabatic vacuum. Thus there is always an inherent ambiguity in the defini- 
tion of the particle concept, and as long as the rate of particle production is 
not negligible, the notion of particle number is not useful. Nevertheless, the 
adiabatic vacuum is the best candidate to represent the vacuum for problems 
that are not asymptotically static because there is a large probability that it 
leaves the high-energy modes empty. This will prove to be crucial, since it 
enables the isolation and regularization of the infinities which appear in the 
expectation values of the currents. A correct way to describe the dynamical 
evolution of the state of a field is to calculate the evolution of local field ob- 
servables, such as {%l)\j^{x)\ip) or (ip\T tiU (x)\tp) , instead of global quantities 
like the particle number. 

The expectation value of the current and adiabatic regularization 

We now turn to the calculation of the current expectation value. As a result 
of spatial homogeneity we have 

(4a q > = (2tt) V(k - q)AL(k) 

(&L k &_ q ) = (27r)V(k-q)AL(k) (2.23) 

(6-ka q > = (2vr)V(k-q)F(k) 

where N+ and AL are the positive and negative charge densities in the initial 
state and F(k) is the correlation pair density or fluctuation term. Gauss' Law 
dictates AL = AL = N, since 

(j°> = e J [dk][JV + (k) - iV_(k)] = 0. (2.24) 
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Using (2.4), (2.9), (2.15), and (2.23), we have 



= e 



(f) 

(k 3 - eA) 



J [dk] ( - 6 1 + 2iV(k) + 2F(k) cos ^2 J n k (t')dt' 



(2.25) 

Here N(k) > and F(k) can be chosen to be real with F(k) > by adjusting 
the overall phase of the mode amplitudes fk(t). The expression in the square 
brackets in (2.25) is positive. To prove this, we note 

< ((4 - 6_ k )(a k - bl k )) = 1 + 2N(k) - 2F(k) , (2.26) 

and hence 

< 2F(k) < 2iV(k) + 1 . (2.27) 

By inspecting the brackets of (2.25), we can identify the first term as the 
vacuum contribution, the second term as roughly the on-shell classical part 
of the current, and the last term as the quantum fluctuation part, which 
oscillates with frequencies larger than 2m. 

Because the densities of particles and correlated pairs at t — are finite, 
only the vacuum term 

e f[dk}{k 3 -eA)^- (2.28) 

contains infinities. In 1+1 dimensions, there is only a normal-ordering infinity 
in the current, which is easily subtracted. In 3+1 dimensions, an infinite 
charge renormalization is needed as well. These divergences can be identified 
and subtracted with adiabatic regularization. 8 

In order to isolate the divergent part of (2.25), we examine the large- 
momentum behavior of Q^. It is convenient to expand 1/f^k m powers of 
1/c^k- (As discussed above, this asymptotic expansion is valid for positive- 
frequency modes.) This is done by successive iteration of the mode equation 
(2.16): Inserting the zeroth-order solution Q^] = uj^ into the right-hand side 
of (2.16) one obtains f2 2 up to second order; inserting this value in the right- 
hand side the fourth order is obtained, and so on. It is not difficult to see 



8 It may be shown [5, 9, 11, 18, 19] that all the divergences obtained in adiabatic 
regularization are the same as those found via point-splitting regularization or via the 
closely related Pauli-Villars [18] and n-wave [4, 5] schemes. 
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that higher-order adiabatic approximations contain terms of higher order in 
l/cuk- We shall see that adiabatic expansion up to second order is sufficient 
to isolate the divergences in the current. Noting that 

^ = - ^ + O(l/o*) , (2.29) 

we have up to second order 

+ 0{l/ul). (2.30) 



Formally we can write the integrand of (2.25) as 

1 eA(k 3 - eA) 



1 + 2N(k) + 2F(k) cos ^2 J* Q k {t')dt' 



1 

(2-31) 

where r is a remainder term which falls off faster than k 4 for k — > oo. By 
substituting (2.31) into (2.25) and replacing k 3 — eA by k 3 in the second term 
of the integral we find 

A = e [ [dk]—^—(k 3 - eA) - e 2 5e 2 A + e [ [dk}(k 3 - eA)r k (t) , (2.32) 
where 

The first term vanishes upon symmetric integration, while the second term is 
the same as the divergent part of the scalar vacuum polarization H(q 2 = 0), 
which is logarithmically divergent in 3+1 dimensions. 

Note that shifting the variable of integration in the first term of (2.32) 
is not trivial, because the integral is linearly divergent when d — 1 and cu- 
bically divergent when d — 3. A correct, gauge invariant way to perform 
the integration and discard this infinity is by cutting off k l at +A + eA 1 , 
which corresponds to fixed integration boundaries for the kinetic momentum 
k — eA. This automatically gives the result of symmetric integration. Alter- 
natively, the infinity can be removed by Pauli-Villars regularization, adding 
one auxiliary massive field to the theory. 
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After removing the first term, we define the renormalized charge and 
gauge field via 

4 = e 2 (l + e 2 5e 2 y 1 = Ze 2 , A R = Z^ 2 A , (2.34) 
and multiply (2.25) by the factor Ze/e R . The renormalized equation is 

A R = e R J [dk] (fc3 ~^ rAr) 1 + 2iV(k) + 2F(k) cos (2 f fl k (t')dt^ 

+e 2 R A R 5e 2 . (2.35) 

The last term formally cancels the divergence of the first. This can be seen 
explicitly by using (2.31) in (2.35), giving 

A R = e R j [dk}(k 3 - e R A R )r k (t) (2.36) 

in terms of the renormalized quantities. Superficially, it seems that the second 
derivative of A appears only on the left-hand side of (2.36), but in fact the 
subsidiary condition (2.31) defining is an intrinsic part of (2.36). 

Since eA = e R A R and since r(k, t) depends on e and A only through the 
combination eA, we will henceforth drop the subscript R. 



Practicalities of numerical solution 

Now that we have renormalized the equations we turn to the practical aspects 
of solving the back-reaction problem. We are interested in solving an initial- 
value problem. The initial conditions for the Maxwell equation are 

A(t = 0) = -Eq , A(t = 0) = . (2.37) 

The initial state of the boson field is the adiabatic vacuum, selected by match- 
ing the exact solutions of (2.10) to their adiabatic approximation, viz., 

n k (t = 0) = cu k (t = 0) , Si k (t = 0) = ufc(t = 0) , (2.38) 

and by setting 

N(k) = F(k) = . (2.39) 

Nonvacuum initial conditions may be specified by adding nonzero particle- 
number and pair-correlation densities to the current expectation value, with- 
out changing the initial conditions (2.38). It is important to note that (2.31), 
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in which r k must fall off faster that k~ 4 at large k, constrains the initial con- 
ditions placed on f2 k . The last two terms in (2.31) are connected at the same 
time by the Maxwell equation (2.36). To check the consistency of the initial 
conditions (2.38), we begin by noting that requiring f2 k = turns (2.31) 
into 

which is a homogeneous equation in A and r k , as is the Maxwell equation 
(2.36). The solution to the system (2.40) and (2.36) is, obviously, A = r k = 0. 
The consistency test comes from differentiating (2.31) and (2.36) with respect 
to time. Upon inserting the initial conditions (2.37) and (2.38), the two 
equations are homogeneous in d 3 A/dt 3 and r k , and possess the consistent 
solution d 3 A/dt 3 = r k = 0. 

Given such a set of consistent initial conditions, one can evolve the set of 
back-reaction equations (2.16), (2.31), and (2.36). By stepping the equations 
forward in time, one arrives at {A, A, f2 k , fi k } at each time. In principle, A 
can be calculated without using (2.36). One can examine the asymptotic be- 
havior of f2(k) at large momentum and extract both r k and A. For instance, 
for large k one can parametrize the left-hand side of (2.31) as 

^"^+4 + 4 + ..., (2.41) 



(for N = F = 0) and directly extract the value of A. The fitted value of A 
depends, however, on the number of terms used to expand the above series. In 
addition, at some momentum, a fit to a finite series in powers of cJ k 1 breaks 
down. Practically one faces another difficulty in solving the semiclassical 
system (2.16), (2.31), and (2.36): A numerical solution of the mode equation 
using a standard scheme, such as Runge-Kutta, does not necessarily agree 
with the correct adiabatic expansion for large momentum, simply because 
the Taylor expansion in the time parameter does not converge for large k. 

Instead, we adopt an iterative scheme making use of (2.36). We take r k 
to be zero at an extremely large momentum as a trial value, so that A is 
extracted from (2.31) automatically. With this value for A we use (2.31) to 
extract r k for each k up to this very large momentum. Then, by substituting 
r k in (2.36) we get a new and slightly different value for A. This procedure 
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may be iterated until convergence of the sequence of (A, r^) is reached. 
Thereafter the next time step is taken. 



2.3 Numerical results in (1+1) dimensions 

It is useful to investigate the set of coupled Klein- Gordon and Maxwell equa- 
tions in 1+1 dimensions, since the momentum integral in (2.36) is only one- 
dimensional and thus the computations are more manageable. The renor- 
malization procedure is somewhat simpler, since the charge renormalization 
is finite [see (2.33)]. Returning to (2.35), we can now solve for A, obtaining 



(2.42) 



The subtraction of l/tVk(t) in the integrand has been inserted by hand. Just 
like the apparent divergence in the unsubtracted integral, the integrated sub- 
traction also vanishes upon symmetric integration. This term is inserted to 
ensure convergence, irrespective of how the numerical integration treats the 
large \k\ regions. Eq. (2.42) is to be coupled to the one-dimensional version 
of (2.16), 

2 



The initial conditions are given by (2.37)-(2.39). 

We show in Fig. 1 the time evolution of the scaled electric field 
E = eE/m 2 and induced current j = ej/m 3 , as functions of r = mt. With 
strong initial electric fields we find that the induced current increases rapidly 
and becomes saturated at a constant value for some time, after which plasma 
oscillations are clearly seen. The fine oscillations in j are not numerical arti- 
facts, since they persist with small time steps (dr = 1CT 4 ) and fine momen- 
tum grids (dk = dk/m = 0.004). 

In order to have some insight into these results, we examine an analo- 
gous back-reaction problem of a classical system of particles and antiparticles 
which interact with a homogeneous electric field. For simplicity, we assume 
that the momentum distributions of the particles and the antiparticles 
are localized at ±p', 

A^Qo, x) = 2n5(p T p')n , (2.44) 
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Figure 1: (a) Time evolution of scaled electric field E and current J, with 
initial value E = 1.0 and coupling e 2 /m 2 = 0.1. Solid line is semiclassical 
scalar QED, and dashed line is Boltzmann-Vlasov model without the stimu- 
lated pair creation correction, Eq. (2.60). (b) Same as (a), but for E = 4.0. 
The dashed curve is as in (a), while the dashed-dotted curve includes the 
correction, Eq. (2.67). 



where n is the density of particles and p is the kinetic momentum. The set 
of coupled equations for this system reads 

A = j = 2env = 2en 



\Jp 2 + m 2 

p = eE = -eA , (2.45) 

where 

p = imv , 1 = {I - v 2 )- 1 / 2 . (2.46) 

The factor 2 in the current accounts for the antiparticles. This system has 
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an integral of motion from energy conservation, 

e = -E 2 + 2n^p 2 + m 2 (2.47) 

and thus 



P 

e = 



k 2 + 2n^p 2 + m 2 , (2.48) 

where e is a constant. The last equation is equivalent to the system (2.45), 
and it reduces the system to the single equation 



^ = (2e 2 e - Ane 2 ^p 2 + m 2 ) ^ . (2.49) 



The solution of this equation is oscillatory with period 

T = 4 / P . (2.50) 

(2e 2 e - 4ne 2 Vp 2 + m 2 ) 

Here p max is defined by the zero of the denominator. In the case of a strong 
initial electric field, the particles and antiparticles are accelerated until they 
approach the velocity of light. At this point the electric current stops in- 
creasing, and the electric field degrades until it changes its direction and 
accelerates the particles in the opposite direction. If the initial field is very 
intense, then most of the time the field is strong enough to keep the particles 
flowing almost at the speed of light. In this case we are able to estimate the 
period T, since the particles are ultrarelativistic. In this limit p 2 ^> m 2 and 
therefore (2.50) yields 



2\ nm + E 2 

T UR = -¥ — , (2.51) 

ne 

where E mSLX is the amplitude of the electric field. The current saturates at 

= 2en\v max \ = 2en . (2.52) 

If the initial electric field is weak, the problem becomes nonrelativistic and 
the set (2.45) is reduced to a simple harmonic oscillator equation with a 
frequency 

2ne 2 m 2tt pm 



m e V 2n 



Up = —i T NR = — X /—. (2.53) 
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Figure 2: Scaled current j at the peak of the first oscillation, as a function 
of the initial scaled electric field E. 

This is the well-known nonrelativistic plasma frequency. It is evident that 
for a fixed particle density n, the plasma frequency decreases as the electric 
field increases. 

The saturation of the first oscillation in our calculation is thus easy to 
understand. At the very beginning of the evolution, particles and antiparti- 
cles are created and accelerated by the electric field. The current j saturates 
as v is driven to the speed of light. Eq. (2.52) gives us a method for estimat- 
ing the number of pairs n created in the first rise of j, allowing us to evade 
the ambiguities in the definition of particle number inherent in the adiabatic 
approach. We display a graph of the peak current in the first oscillation as 
a function of the initial field in Fig. 2. 

In subsequent oscillations, n is larger and the peak value of E is weaker 
because of particle production, so the particles remain nonrelativistic even at 
the peak of the current. At late times, when the envelope of the electric field 
approaches a constant, that is, when further particle production is negligible, 
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one can estimate the number of particles from (2.53). 

We can also evaluate the density of particles by calculating the number 
operator at some late time, with the expectation that at late times the num- 
ber density is almost conserved. We define the particle number density by 
expanding the exact bosonic field operator $ in terms of the zeroth-order 
mode functions 



fl°\t) = 



1 



V^k(*) 



cxp 



-% y* uj k (t')dt' 



(2.54) 



Then the corresponding creation and annihilation operators, and b k \ 
become time-dependent, and the vacuum expectation value of the number 
operator 

n(k;t)^(0\a^(t)af\t)\0) (2.55) 

may be computed by a Bogolyubov transformation from the basis functions 
(2.15) to the basis functions (2.54). (See Appendix A for details.) In the 
limit of large times we define 



lim a ( k °\t) =a° k ut . 

t — >oo 

The expectation value of the number operator for large times reads 



(2.56) 



n(k) = (0\a°^ar\0) 



= lim 



a k 
1 



t^oo AVl k uj k 



(Q k - uj k ) 2 + - ((fifc/Qfc) - {uJk/uJk 



(2.57) 



We plot this quantity in Fig. 3. 

The rapid oscillations in n(k) are related to rapid oscillations in the in- 
tegrand in (2.42), which develop with time. From (2.17) we obtain 



n k (t) 



cik\ 2 + \c 2 k\ 2 + c* lk c 2 k exp (2i J W k (t')dt' 



W k (t) 

+ c lk c* 2k exp (-2i J W k {t')dt' 



(2.58) 



Since W k (t) ~ u k (t) in the absence of interaction [see the discussion after 
Eq. (2.19)], it is evident from (2.58) that for a fixed time, fl^ 1 oscillates as a 
function of k. At a fixed time t, this expression oscillates in momentum space 



21 




Figure 3: (a) Momentum distribution of produced pairs, for the evolution 
shown in Fig. 1(a), at time r = 550. The abscissa is the scaled kinetic 
momentum p = k — A, with k = k/m. (b) Data of (a) after smoothing 
(histogram), compared with Boltzmann-Vlasov model. The smooth curve 
is for the model with Bose enhancement, the broken curve for the model 
without. 
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with a frequency on the order of t. As time goes by, then, the frequency grows, 
and hence the integrands in the current (2.42) and the particle distribution 
(2.57) also oscillate. This kind of oscillation appears even if the initial pair- 
correlation density (related to zitterbewegung [43]) is zero. 

When the interaction is switched off, Co^ = but Cl^ ^ 0. It is clear that 
in the absence of interaction the expectation value of the number operator 
should be a constant. One finds that the particle density in (2.57) is indeed 
time-independent for tOk = by calculating h(k) using (2.16). 



2.4 Phenomenological kinetic approach 

The simple classical picture presented in the preceding subsection can ex- 
plain gross features of the time evolution of our system, such as the early 
saturation of the current and, for that matter, the very existence of plasma 
oscillations. Particle creation is not included, nor is the distribution of the 
created particles in momentum. The inclusion of these features requires a 
more sophisticated model, and a relativistic kinetic equation suggests itself. 
The original Boltzmann-Vlasov equation, however, conserves particle num- 
ber, and thus one must add an explicit source term. 

Many calculations in the realm of nucleus-nucleus collisions have been 
based on such a phenomenological model [39-42]. In the case of a homoge- 
neous electric field in 1+1 dimensions, the relativistic kinetic equation is 

df | cE df _ dN 
dt dp dtdxdp 

where f(p,t) is the (x-independent) classical phase-space distribution, ex- 
pressed as a function of the kinetic momentum p, and the right-hand side is 
the boson pair-production rate. For the latter, one assumes the applicability 
of Schwinger's formula, 



\eE(t)\\og 



dt dx dp 



I Tim 

1 + exp 



2 \ 



\eE(t)\J 



S(p) , (2.60) 



despite the fact that here the electric field is not constant in time. The 
factor 5(p) indicates the usual WKB assumption that the particles are pro- 
duced at rest [30]. Initially, f(p,0) = 0. Eq. (2.59) may be solved using the 
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characteristics dp/dt = eE, giving 
f(p,t) = 



\eE(t') 



dt'\eE(t')\ log 1 + exp 
x S(p - eA(t') + eA(t)) . 
The 5-function allows us to perform the integration, whence 



/(p,0 = E 1o s 



1 + exp 



Km 



\eE(U 



(2.61) 



(2.62) 



where the Vs fulfill p + eA(t) - eA(U) = and U < t. 

The kinetic equation is coupled to the Maxwell equation, 



£A 
dt 2 



Jtotal Jcond Jpol i 



where the conduction current is 



j cond = 2e J T^^-fiP't) > 



(2.63) 



(2.64) 



with e p = y/p 2 + m 2 , and the polarization current is [42] 



Jpol 



dp 



dN 



E J 2vr p dtdxdp 



(2.65) 



[The factors of 2 in (2.64) and (2.65) account for the contributions of the 
antiparticles.] Inserting (2.61) into (2.63) reduces the system to a single 
equation, 



d 2 A 
dr 2 



-Km- 



[ T dr' 
Jo 

x log 



Mr') ~ Mr) 



At?) - A(r 
1 + exp 



+ 1 



E(r') 



IT 



+ 



nm- 



■ sign (.E(r)) log 



E(t>) 
1 + exp 



K 



E(t) 



(2.66) 
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in terms of the dimensionless variables A = eA/m, E, and r. 

Equation (2.59), widely used in the literature, omits a statistical factor 
which should be present. Roughly speaking, the source term should contain 
the Bose-Einstein factor (1 + f(p,t))(l + f(—p,t)), where / is the antipar- 
ticle density. The CP symmetry of the problem (and of the chosen initial 
conditions) implies that f{—p,t) = f(p,t), and thus the statistical factor is 
(1 + f(p,t)) 2 . Detailed balance, however, dictates that there should be a 
loss term as well, due to particle annihilation, containing the same matrix 
element but the statistical factor f(p,t)f(—p,t) = (f(p,t)) 2 . The difference 
of the gain and loss terms is thus proportional 9 to (1 + f(p, t)) 2 — (f(p, t)) 2 = 
1 + 2f(p,t). We thus replace (2.60) with 



dN 

:i + 2/(p,*))|e£(i)|log 



dt dx dp 



I Tim 

1 + exp 



2 \ 



\eE(t)\ 



S{p) ■ (2.67) 



The modifications to (2.61)-(2.62) and (2.66) are obvious. In (2.62), note 
carefully the strict inequality ti < t. 



2.5 Comparison in 1+1 dimensions 

Let us compare the results of the phenomenological kinetic model with those 
of our semiclassical analysis. The time evolution of E and j in the kinetic 
theory is shown in the dashed and dotted curves of Fig. 1, where we see that 
there is good quantitative agreement between the results obtained with the 
two very different methods, as long as Bose-Einstein enhancement is properly 
taken into account. 10 Still, the oscillations are faster and the electric fields 
decay more rapidly in the quantum calculation than in the Boltzmann-Vlasov 
model. 

The distribution function f(p,t), measured after several plasma oscilla- 
tions, may be compared with the n{k) of the quantum theory after the latter 
is smoothed, as shown in Fig. 3. Naturally, the curves have different normal- 
izations and a relative displacement due to the slightly different values of j 
and A. The kinetic theory is quite successful at reproducing the (smoothed) 

9 We thank J.-P. Blaizot for this simple argument. For a full derivation of the statistical 
factor see [51]. 

10 The enhancement is not significant in the case of Fig. 1(a), where the field is compar- 
atively weak and particle production is slow. 
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quantum result. One can thus use the kinetic-theory model to explain various 
features of the particle distribution, such as the sharp edges and tails. 

Let us investigate the formal solution given in (2.62). The sum in (2.62) 
consists of terms for which the following condition is fulfilled: 

p + eA(t) - eA(ti) = 0, t t < t. (2.68) 

This is represented graphically in Fig. 4. Eq. (2.68) shows that the maximal 
and minimal values of A at all times up to t determine the highest and lowest 
momenta in the momentum distribution function, since for p > eA ma , x — eA(t) 
and for p < eA m { n — eA(t) there is no solution for t{. Now consider Fig. 5. 
In the case where E t= q = 4 the A field is always negative, never crossing 
zero after r = 0. By examining Fig. 5(b) we can explain the existence of the 
plateau on the right-hand side of the distribution function n{k): For momenta 
in the range eA^l x — eA(t) < p < eA(0) — eA(t), where A^l x is defined in 
Fig. 4, the only contribution to the sum (2.62) results from ti near zero, and 
the electric field at these early times is changing slowly from its initial value. 
The momentum distribution drops suddenly to zero at p — eA(0) — eA(t), 
where the only solution to (2.68) disappears and there are thus no terms in 
the sum (2.62). 

Consider now a generic point p in the momentum distribution at time t. 
There are several tj's which fulfill condition (2.68) for p = p. If we increase p, 
there are branch points p* at which the number of tj's is reduced by two, and 
therefore there are two fewer terms in the sum (2.62). For these values of p, 
the dashed line in Fig. 4(b) just touches a local extremum of A. As p — > p*, 
two of the times ti approach t*, the time at the extremum of A. Note, now, 
that E(t*) = —A(t*) = 0. Terms that have E(ti) = do not contribute to 
the sum (2.62) and therefore at these branch points in p the distribution is 
continuous, though its derivative jumps. 

There are, however, two special cases. The point p which corresponds 
to ti — is one — here, E(t{) is the initial value of the electric field, which 
is nonzero. The momentum distribution is thus always discontinuous at 
p = eA(0) — eA(t). The right edge of the distribution in Fig. 5(b) is a special 
case of this; in Fig. 5(a), the discontinuity is at p ~ 8, in the interior of the 
distribution. 

The other special case is p = 0. Here we have the disappearance of the 
solution of (2.68) with ti = t, at which, again, E ^ 0. Thus there is always 
a discontinuity at p — 0. 
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Figure 4: (a) Schematic plot of A(t), showing graphical construction which 
gives f(p, t) from (2.62). The boxes represent the set of ti corresponding to a 
given value of p, which is the length of the vertical arrow at t. We define 
to be the value of A(t) at the highest stationary point, (b) Construction of 
f(p,t) when p = p*, a branch point where terms in (2.62) disappear. 
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Figure 5: (a) Time evolution of scaled gauge field A, electric field E, and 
current j for the Boltzmann-Vlasov model with the stimulated pair creation 
correction. The initial field is E = 1.0 and the coupling is e 2 /m 2 = 0.1. Also 
shown is the momentum distribution / of produced pairs at time r = 550. 
(b) Same as (a), but for initial E = 4.0, with / measured at r = 80. 
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In the case of -E'(O) = 1, shown in Fig. 5(a), A^ x > A(0), and the 
distribution function does not vanish at the point p = eA(ti = 0) — eA(t). 
Rather it decays smoothly to zero as p — > eA^l x — eA(t). 

We offer a final observation regarding the Bose factor 1 + 2/. If the 
electric field were never to change sign, the factor would have no effect, since 
new particles are only created at p — 0, while the old particles have been 
accelerated by the field. Only if E changes sign do some of the existing 
particles return to p = so that f(p = 0) is nonzero. Thus the Bose factor 
is not relevant for the usual case of the Schwinger mechanism in which the 
electric field is fixed in time. 

2.6 Numerical results in 3+1 dimensions 
Relation to the flux-tube model 

So far we have considered calculations in 1 + 1 dimensions. Turning to 3+1 
dimensions, let us begin by fixing the parameters and initial conditions so as 
to correspond to the physical problem of particle creation in the flux tube 
model [30, 39, 40]. The string tension a, interpreted as the energy per unit 
length stored in the flux tube between quark and antiquark, is given by 

a = X -E 2 A , (2.69) 

where A is the cross-sectional area of the flux tube and E is the strength of 
the longitudinal electric field. Applying Gauss' Law to the quark at either 
end of the flux tube, one obtains 

EA = e, (2.70) 

where e is the effective charge of a quark (sweeping non-Abelian features un- 
der the rug). The Regge slope gives for the string tension [30] a ~ 0.2 GeV 2 . 
Taking 2.5 GeV -1 for the radius of the tube, and using constituent quark 
masses of 0.35 GeV, one finds by using (2.69) and (2.70) that the dimension- 
less field strength in this elementary flux tube is E = eE/m 2 ~ 3, and the 
charge is e ~ 2.5. The strength of the chromoelectric field created in high- 
energy nuclear collisions should be stronger than this, since the sources will 
be in higher representations of the color group [55, 56]; in our calculations, 
therefore, we choose initial fields in the range E =6-10, with e 2 =4-10. 
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Numerical results 

We solve the system (2.16), (2.31), and (2.36) by using Runge-Kutta along 
with the iterative scheme described at the end of subsection 2.2. Previous 
attempts to apply Runge-Kutta to the mode equation have encountered in- 
stabilities [58], but we find that the combination of Runge-Kutta with the 
iterative scheme proves to be stable. Simpson's rule is used to perform the 
integration in (2.36) in cylindrical coordinates. 

In Fig. 6 we display the dimensionless gauge field A, electric field E , and 
current j as functions of the dimensionless time r. We have set -E'(O) = 6 
and e 2 = 10, and measure energy in units of m. To get these high-precision 
results, we used a momentum matrix grid of 200 x 1000 for transverse and 
longitudinal momenta, with a time step dr = 0.0005 and momentum intervals 
dk\\ = dk± = 2n/25. One can see the saturation of the current in the first 
oscillation, as in 1+1 dimensions. To go to large times, we used smaller 
momentum grids and larger time steps. For the results shown in Figs. 7 
and 8 we used momentum grids of 100 x 400 and 200 x 200, respectively, 
and a time step of dr = 0.005, with momentum intervals as in the case shown 
in Fig. 6. Here we can see plasma oscillations. 

In order to compare the semiclassical calculation to the phenomenological 
model, we simply change m 2 to m\ = m 2 + k\ in (2.59)-(2.62) and change 
the phase space in (2.64)-(2.65) to rf 3 p/(27r) 3 . The result in 3 + 1 dimensions 
corresponding to (2.66) is 



including Bose enhancement; the tilde represents dimensionless variables. In 




x 



A{t>) - A{t) 





(2.71) 
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Figure 6: Time evolution of scaled gauge field A, electric field E, and current 
j in 3 + 1 dimensions, with initial value E = 6.0 and coupling e 2 = 10. 
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Figure 7: Time evolution of scaled electric field E and current J in 3 + f 
dimensions, with initial value E = 7.0 and coupling e 2 = 4. Solid line is 
semiclassical scalar QED, dashed line is the Boltzmann-Vlasov model with 
the stimulated pair-creation correction. 
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Figs. 7 and 8 it is seen that the phenomenological Boltzmann-Vlasov re- 
sults agree fairly well with those of semiclassical QED. Still, the quantitative 
agreement is not as good as in 1+1 dimensions; furthermore, the phenomeno- 
logical model misses entirely the rapid oscillations superposed on the plasma 
oscillation in J(r). 

3 Fermion pair production 

3.1 Introduction 

We present here the extension of the semiclassical formalism to spin-| fields 
[51]. Again, we are limited to the case of a spatially homogeneous, classical 
electric field. In the preceding section we showed that the results of our 
calculations in semiclassical scalar electrodynamics are very similar to those 
obtained from the phenomenological relativistic Boltzmann-Vlasov equation. 
It is of interest to see if the same results obtain for fermions which, unlike 
bosons, possess no classical limit, and for which Pauli blocking, in place of 
Bose enhancement, implies quite different consequences. 

To apply adiabatic regularization to the fermion problem, we shall again 
need to express the Fourier amplitudes of the field operators in a WKB-like 
form. This will allow us to isolate the ultraviolet divergences through an 
adiabatic expansion and to eliminate them through renormalization. The 
WKB parametrization will turn out to be a bit more complex than that for 
the bosons. We shall further see that in order to make the current zero in the 
initial adiabatic vacuum, we shall need to express physical expectation values 
by averaging over two different complete sets of linearly independent solutions 
of the Dirac equation. Our initial state is thus a mixed state. In Subsection 

3.2 we derive the coupled equations for the fields in the semiclassical limit of 
QED, and discuss the adiabatic regularization procedure. In Subsection 3.3 
we present our numerical results in 1+1 dimensions and compare them with 
the phenomenological Boltzmann-Vlasov model. 
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3.2 Fermion QED in the semiclassical mean field ap- 
proximation 

Equations of motion 

The lagrangian density for electrodynamics is 

C = m^{d^ + ieA^ - - ^F^F^, (3.1) 

where the metric convention is taken as (H ). We work explicitly in 3+1 

dimensions, leaving the case of 1+1 dimensions for later. For the 7-matrices 
we use the convention of Bjorken and Drell [59], 

where I is the identity matrix and a 1 are the Pauli matrices. 

Again we quantize only the matter field, while the electromagnetic field 
is treated classically. The coupled field equations read 

(i^d^ - e^A^ - m)V(x) = (3.3) 

and 

d,F^ = (f) = ^,^]) , (3-4) 

where the expectation value is with respect to the initial state of the spinor 
field. The commutator in the electric current guarantees a zero expectation 
value for any charge-conjugation eigenstate. Expressing the solution of the 
Dirac equation as 

q?( x ) = (i-fd,* - e 7 M M + m)Q(x) , (3.5) 

and inserting (3.5) into (3.3), it follows that $ satisfies the quadratic Dirac 
equation [60] 



(id ll -eArf--a^F^-m 7 



= , (3.6) 



where $ is a four-component spinor. Here we consider the case where the 
electric field is spatially homogeneous so that the field strength F^ v depends 
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only on time. Owing to homogeneity, the Maxwell equations (3.4) allow 
only configurations where = 0. We again take the electric field to be in 
the direction of the z-axis and choose a gauge such that only A = A 3 (t) is 
nonvanishing. Then the second-order Dirac equation becomes 



□ + e 2 A\t) + 2iA(t)d 3 - ie<9 A(t)7°7 3 



m 



®(x) 



0. 



(3.7) 



Consider first the c-number Dirac equation corresponding to (3.7). Spa- 
tial homogeneity implies that there exist solutions of the form 



<Pks(x) = e lk - x / ks (t)x s , 



where 



Xi 



X3 



( 1 \ 


1 

V o / 

/ i 
o 

-l 

V o J 



X2 



X4 



/ \ 
1 



\ -1 J 

(°\ 
1 



V i / 



(3.8) 



(3.9) 



(3.10) 



These spinors are chosen to be eigenvectors of a 3 = 7°7 3 in the representation 
(3.2) for the 7-matrices — the eigenvalues of a 3 are X s — 1 for s — 1,2 and 
A s = — 1 for s = 3, 4. The spinors satisfy the normalization and completeness 
conditions 



T,(xtr(Xs)a = 25 rs , 



a=l 
4 



Eixinxr), = 2S« 



p ■ 



(3.11) 



r=l 



Substituting (3.8) into (3.7), it follows that the mode functions fk s (t) satisfy 



d 2 hs(t) 
dt 2 



+ 



dA 
iKe— 
dt 



hs(t) = , 



(3.12) 



with 



u£(t) =pl + k 2 ± + m 2 , k 2 ± = k 2 + kl pi^tf-eA^t) 



(3.13) 
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Equations (3.12) are second-order differential equations, and therefore for 
each s there are two independent solutions. Let f£ s and f^ s be the two inde- 
pendent solutions of (3.12), which become positive- and negative- frequency 
solutions in the absence of the electric field. Clearly at the moment we have 
eight different solutions for the second-order equation (3.6), namely, ff s for 
s = 1,2,3,4. However, the Dirac equation (3.3) has only four independent 
solutions. If we restrict ourselves to solutions that belong to the set s = 1,2 
or to the set s = 3, 4, we shall see that from each set one can construct a 
linearly independent set of solutions of the Dirac equation. 

The form introduced in (3.8) allows us to write ip as 

^ ks (x) = Dfasix) = (i^do + - e 7 3 A 3 + m)(j) ks (x) . (3.14) 

Explicitly, the two sets of independent solutions of the Dirac equation may 
be taken to be 

rl>£ = e ik - x Df± X ., s = l,2, (3.15) 

and 

tf£ = J^DftXs , * = 3,4. (3.16) 
Using Eqs. (3.14)-(3.16) we find, for a given k, 

-*\sP3 [ffff - if] } Srs , (3.17) 

where either r, s — 1, 2 or r, s — 3, 4. An exactly analogous formula may be 
derived for ip^t/jj. By differentiating these expressions with respect to time 
and by using Eq. (3.12), it can be readily verified that these inner products are 
time-independent. As t — > — oo there is no interaction between the fermion 
field and the electromagnetic field, and we can choose two independent plane- 
wave solutions for Eq. (3.12), 

lim f£ = c a e***\ (3.18) 

t — >— oo 

where c s are constants. Insertion of these free solutions in the relation for 
ip^ipj yields immediately 

^V? = . (3.19) 
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Since this result is time-independent, it is valid at any time, and each set ipf, 
with s = 1, 2 or s = 3, 4, is a complete set of linearly-independent solutions 
of the Dirac equation. Note that these complete systems are not identical, 
and orthonormality conditions holds for each set separately. In principle, we 
need only one of these sets in order to expand the field operator \1/ in terms of 
single-particle solutions. In order to ensure that, with our initial conditions, 
the Dirac current vanishes at t — 0, it is advantageous to use both sets in our 
calculations. 

We now construct the quantized spinor field operator in the form 



*(*) = f[dk] E + ^HOifc] 

J s=l,2 

= / [dk] J2 [b s (k)^ s + dt(-k)^,] , (3.20) 

J s=3,4 

where the two lines show the field expressed in terms of the two alternative 
bases. Here [dk] = c? 3 k/(27r) 3 . The fermion fields obey canonical anticom- 
mutation relations, {^ a (t,x), ^g(t, y)} = <5 3 (x — y)5 a g. The creation and 
annihilation operators of each set, r, s — 1, 2 or r, s — 3, 4, obey the standard 
anticommutation relations, 

{b r (k), &t(q)} = K(k), rft(q)} = (27r) 3 5 3 (k - q)5 rs , (3.21) 

if we impose the condition 

^Vf = $rs , (3.22) 

to normalize the mode functions. 

We now calculate the expectation value of the electric current. For the 
sake of simplicity we choose the initial state to be the vacuum annihilated 
by b r (k) and d r (k). Using the anticommutation relations (3.21) we find 

(0|J 3 |0) = ^<0|[^,7 3 *]|0> 

= | / [dk] E {-^ f 7°7 3 ^ + ^slW^s} ■ (3-23) 

Z J s=l,2 

Alternatively, 

(0|J 3 |0) = I I [dk] E {-^W^ + SVfe } • (3.24) 

Z J s=3,4 
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Averaging the two expressions, 



(0|j 3 |0) = - ( [dk] £ {-^ f 7°7 3 ^ + Wife } • (3.25) 

4 J s=l 

This form will be useful when we turn to the adiabatic regular izat ion of the 
current. The other components of the current are zero since the electric field 
is in the z-direction. Using (3.14) we find 



0„,3„/,± 
ks 



A, 



-ij°d + 7 X k ± - (k 3 - eA 3 )7 3 + m 



(3.26) 



and thus (3.17) and (3.22) give 



2 i 2 

1 m 



P3 2 )f£ft - tiftis + ***(f£fc - id Jus)} 



= A, 



(3.27) 



where the index s is not summed over. Inserting (3.27) into (3.25) yields 



(0|J 3 |0) = e£ / [dk](fci + m 2 )A s (\f k f - |/+| 2 ) . 

s=l J 

From (3.17) and (3.22) it can be shown [60] that 



2(A:i + m 2 )(|/+| 2 + |/ 1 
Eq. (3.29) then gives the current as 



ksl 



1. 



(0\f\0) = -2ej:f[dk](kl + m 2 )X s \f ] 

s=l J 



+ 1 2 

ksl • 



(3.28) 



(3.29) 



(3.30) 



Averaging the current expectation value over the two bases as in (3.25) means 
that we prepare the initial system in a mixed state. 



Adiabatic regularization 

As in the boson case of the previous section, the difficulty in solving the 
coupled semiclassical equations (3.3) and (3.4) originates from the fact that 
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the expectation value of the current (3.30) diverges in the interacting theory. 
This infinity can, as before, be removed by charge renormalization. In order 
to isolate the ultraviolet behavior of the current integrand by adiabatic ex- 
pansion, we need to express the mode equations (3.12) in a WKB-like form. 
The generic problem is to find a suitable parametrization for the solution of 
the differential equation il(t) + e(t)u(t) = 0, where in the present case e is 
the complex quantity in square brackets in (3.12). Such a parametrization is 
found in [57], namely, 

/i(t) = ^ vm exp K (- iSi "« - a -2?£m) *'} • (3 - 31) 

where are normalization constants and f2ks is a real generalized frequen- 
cy. 11 By substituting (3.31) into (3.12) we obtain the WKB-like equation for 



n^(0 = - — + - ^ + o£ + ( eA V - \ s — + \ a e -^. (3.32) 



As in the boson case, the equation for Q is a second-order nonlinear differ- 
ential equation. 

This equation enables us to study the large-momentum behavior of the 
solutions. As above, an adiabatic expansion of (3.32) to second order is 
needed to identify the divergences in the current (3.30). Noting that 

„ k = z eAp 1 + 0(1/cJk)) (3>33) 

we have up to second order [i.e., iterating (3.32) once] 

fi ks = LJ k -eA (\ s iu k -p 3 ) /Auol + 0(l/c4). (3.34) 
Using the ansatz (3.31) the current reads 



(0|j 3 |0) = -2e£ [[dk](kl + m 2 )\ s 

s=l J 



2^ks 1 J ^ks(t') 



(3.35) 



11 The second solution f^ s for the mode equation can be found by using its Wronskian. 
When choosing the form (3.31) for f£ s , the second solution does not have a simple form, 
and for this reason we expressed the current (3.30) in terms of the positive-frequency 
solutions only. 
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Equations (3.17) and (3.22) determine the normalization constants, and it 
turns out that the expression in square brackets is 



ks 



ks 



exp <^ —X s 



eA(t') 



■dt' 



(4 + 



-X.eA 



n 



2Q 



ks 



2Q 



k > + n 2 



ks , 



ks 



- \ s 2p 3 tt ks 



= r.(k) 



With the identity 

we obtain 
(0|J 3 |0) 



k\ + m 2 = [u k + p 3 ] [u k - p 3 ], 



(3.36) 
(3.37) 



e J [dk] | 



^k ~P3 



^k + £>3 



+ ok- 4 )}. 

At large momentum, we approximate 

1 e'A 



{u k + eA(u k + p 3 )/4ul to k - eA{u k - p 3 ) / Auj^ 

(3.38) 



l±eA[cu k ±p 3 }/Acui 



~ 1 =p 



4^ 



[w k ± Ps] 



(3.39) 



After we perform the angular integrations and drop terms that are odd func- 
tions of p 3 , the Maxwell equation becomes 



A = (0\f\0) = -e 2 Aj[dk]^ 
= —e 2 A 5e 2 + (finite part), 



2u| - ||) + ( finite P &rt ) 



(3.40) 



where 



[dk] 



dk 



2ui 2ulJ 47r 2 7o L(P + m 2 )i 3(fc 2 + m 2 



(3.41) 

The current in (3.40) diverges logarithmically, with the same divergence 
as the vacuum polarization Il(g 2 = 0). Renormalizing as usual, we define 



e 2 {l + e 2 5e 2 )' l = Ze 2 , A R = Z' 1 / 2 A, 



(3.42) 
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so that eA = crAr. We can also write Z — (1 — e R 5e 2 ). Multiplying (3.35) 
by Ze/eR we obtain 

AR-e 2 R A R 5e 2 = e R j2 f[dk](k 2 ± + m 2 )(-\ s )T s (k) . (3.43) 

s=l J 

Using (3.41) and rearranging, we have 



A R = e R J [dk] 



(k 2 ± + m 2 )J2(-K)r s (k)+e 2 R A R 



(3.44) 



where the right-hand side is finite. We drop the R subscripts now. 

As for the boson problem, we can define a remainder which is the dif- 
ference between the integrand in (3.35) and its adiabatic approximation. 
Examining (3.38), (3.39), and (3.43), we can write 

(*± + m 2 ) £(-A s )r s (k) = - J* (o£ - n 2 3 ) + i2 k (t) . (3.45) 

s=l ^ k 

At large momentum this minimal adiabatic approximation matches the exact 
integrand up to terms that fall off as 0(l/u%.), so the remainder Rk(t) falls 
off faster. Upon substituting (3.45) into (3.44) and using (3.38)-(3.40), the 
finite Maxwell equation takes the form 



J[dk]R k (t). (3.46) 



Again, the subsidiary condition (3.45) defining i? k is an intrinsic part of 
(3.46). 

As in the boson case, the initial conditions we impose are 

A(t = 0) = -E , A(t = 0) = , (3.47) 

and 

^ ks (t = 0) = u k (t = 0) , (l ks (t = 0) = u k (t = 0) , (3.48) 

which specifies the adiabatic vacuum. Nonvacuum initial conditions may 
be handled in a manner analogous to the bosonic case by adding nonzero 
particle number densities to the current expectation value. 
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Here, too, the choice of initial conditions is constrained by demanding 
the consistency of the adiabatic expansion (3.45) with the Maxwell equation 
(3.46). By substituting (3.48) into (3.36) and (3.45) we find that A(0) = 0, 
but 



is not zero. It turns out, however, that the integration over k in (3.46) gives 
zero by charge conjugation symmetry 12 p 3 — > — p 3 . Thus the initial conditions 
(3.47) and (3.48) are consistent with the requirements of renormalization. 
We emphasize that choosing i?k(0) = is not consistent with the initial 
conditions (3.48). 

Given such a set of consistent initial conditions, we can solve the back- 
reaction equations (3.32), (3.45), and (3.46) exactly as we did in the scalar 
case. To do this we take i?k to be zero at an extremely large momentum 
as a trial value, so that A can be extracted from (3.45). Now having A, we 
use (3.45) to extract i?k for each k up to this very large momentum. Then, 
substituting in (3.46) we get a new corrected value for A. This procedure 
may be iterated until convergence for A and i?k is reached, after which one 
may proceed to the next time step. 

3.3 Calculation in 1+1 dimensions 

Dynamical equations 

For the fermion problem, we present numerical results only for the (1+1)- 
dimensional case. Let us indicate the modifications needed in the equations. 
The 7-matrices are given by 



12 Thc summation in (3.45) over s, along with the A s factor, is essential for making the 
integrand odd in p 3 . It is here that our adoption of a mixed initial state proves its value. 
Choosing only the s = 1,2 basis, for example, would have made it impossible to have a 
zero initial current. 




(3.49) 




(3.50) 
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and 7 1 plays the role that 7 s played in the (3+l)-dimensional case. In A = 
gauge, we define A 1 = A(t), and the second-order Dirac equation is 



□ + e 2 A\t) + 2iA(t)di - ieflb^WYV 



m 



0. 



(3.51) 



The Dirac equation in two dimensions has two independent solutions: Either 



fkiXi , fkiXi 



or 



fk2X2 , fk2X2 

may be taken as the basis set of independent solutions. Here 7°7 1 x* 
with Ai = 1 and A 2 = — 1, and the spinors Xs are given by 



We define 



Xi 



2 2 
p + m , 



X2 



1 
-1 



p 



k-eA\t) . 



The Maxwell equation is 

A = (01/10) 



dk 



2e ^ 2 {\m 2 -\f k \n- 



2tt 



(3.52) 
(3.53) 

= ^sXs 

(3.54) 
(3.55) 
(3.56) 



In two dimensions there is no spin, and thus there are half as many terms as 
in (3.30) after summation. 

Renormalization of (3.56) can be done in the same way as in subsection 
3.3. In two-dimensional QED the charge renormalization is finite and we find 
that 5e 2 = (67rm 2 ) -1 . Therefore, in the renormalized Maxwell equation, A 
can be isolated [in contrast to (3.44)], and we obtain 



en 




where 



r s (k) = 



e 2 R 5e 2 



2 , 1 ^ s zrAr 



m 2 £ (_A s r s (A;)) + 2^ 

.=1,2 W k 



2tt 



ks 



2VL 



+ tt 2 ks - \ s 2ptt 



ks 



ks , 



(3.57) 



(3.58) 



The last term in the braces in (3.57) does not contribute to the integral 
but was included for numerical purposes. The set of equations (3.32) and 
(3.57) with the initial conditions (3.47) and (3.49) defines the numerical 
back-reaction problem. 
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Numerical results 

We show in Figs. 9 and 10 the time evolution of the scaled electric field 
E = eE/m 2 and the induced current J = ej/m 3 as functions of r = mt. 
Here we used a time step of dr = 5 x 10~ 4 and a momentum grid with 
dk = dk/m = 0.003. The results are similar to those for bosons, with 
initial particle creation followed by plasma oscillations. In this case the 
current saturates twice in the early stages of evolution, each time because 
the particle velocity approaches c, with more particles present (and hence a 
larger current) in the second instance. 

The amplitude of the electric-field oscillations decreases substantially only 
in the first few oscillations and remains almost constant at later times. This 
means that essentially all of the pair production happens in the first oscil- 
lations. One sees that in successive oscillations n is larger and E is weaker, 
and therefore the frequency of oscillations increases. As we have noted, the 
frequency of relativistic plasma oscillations depends not only on n but also 
on the amplitude of E: The weaker the field the higher the frequency, in con- 
trast to the nonrelativistic case where the plasma frequency does not depend 
on the amplitude. 

We display in Figs. 11(a) and 12(a) the momentum distribution of the 
produced particles, defined by 13 



^ 3 = 1,2 

In accordance with the Pauli principle, the occupation number does not 
exceed one. 

Phenomenological transport equations 

The fermionic counterpart of (2.67) is 



13 Cf. (2.55); n here is an average over spin states. For the explicit form of n(k;t) for 
large times, analogous to (2.57), see [51]. 




(3.59) 





dN 



(3.60) 
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Figure 9: Fermion pair creation in 1+1 dimensions: Time evolution of the 
scaled electric field E and current j, with initial value E = 1.0 and coupling 
e 2 /m 2 = 0.1. Solid line is semiclassical QED, and dashed line is Boltzmann- 
Vlasov model. 
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Figure 10: As in Fig. 9, but for initial field E = 4.0. 
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Figure 11: (a) Momentum distribution of produced pairs, for the evolution 
shown in Fig. 9, at time r = 600. The abscissa is the scaled kinetic momentum 
p = k — A, with k = k/m. (b) Data of (a) after binning (histogram), 
compared with Boltzmann-Vlasov model (curve). 
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Figure 12: The same as Fig. 11, but for the evolution shown in Fig. 10 (i.e., 
initial field E = 4.0), at time r = 200. 
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Recall that p is the kinetic momentum. The right-hand side is the Schwinger 
expression for the fermion pair-production rate in 1+1 dimensions, multiplied 
by the statistical factor (1 — 2/) which represents Pauli blocking and detailed 
balance (see the brief discussion above Eq. (2.67)). 

Note that the singularity on the right-hand side of (3.60) for m = is 
spurious. The Schwinger formula is based on a tunneling picture [30] that 
incorporates conservation of energy, but ignores other conservation rules that 
would follow from specific features of the system in question. In the present 
case, where dynamical photons are absent, chirality conservation actually 
forbids pair production from the homogeneous electric field for massless par- 
ticles in one spatial dimension, so that the use of the Schwinger formula is 
unjustified in this limit. In the field theory, however, there is no such limita- 
tion. The m — * limit is nothing other than the Schwinger model [28, 61], 
and the chiral anomaly gives a finite rate [28, 62] for pair production. 

Eq. (3.60) is solved in the same way as (2.67). The time evolutions of the 
scaled field strength E and current j are shown in the dashed curves of Figs. 9 
and 10. The quantitative agreement between the kinetic theory and the 
quantum theory is striking. The amplitude of the electric field approaches a 
limiting value after a few oscillations, meaning that thereafter the production 
of particles is negligible. In the boson case an analogous effect is seen, but 
it sets in somewhat later than for fermions. The constant amplitude reflects 
the absence of pair creation from virtual photons and the exponentially small 
spontaneous pair creation rate at this stage of the evolution. The fact that 
the electric field reaches its limiting value more quickly for fermions than 
for bosons may be due to the difficulty of producing more fermions once the 
low-momentum states have been occupied. 

The comparison between f(p,t) at a late time and the smoothed n(k) of 
the quantum theory is shown in Figs. 11(b) and 12(b). 14 The importance 
of Pauli blocking is obvious. Without the statistical factor in (3.60) the 
occupation number exceeds one if the initial electric field is strong enough 
[51]. The Schwinger source term in (3.60) is very large for a strong electric 
field; it is the (1 — 2/) term that prevents the violation of the Pauli principle. 

As for the boson case, our comparison provides a justification for using 
the kinetic theory in studying the quark-gluon plasma. Past treatments have 
not, however, included the Pauli-blocking term which is crucial for strong 

14 Again, the curves have a relative displacement due to the slightly different value of A. 
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fields. 



4 Concluding remarks 

We have solved semiclassical QED for strong fields using the adiabatic reg- 
ularization method for spatially-homogeneous electric fields. We determine 
initial values of the fields and their derivatives which are consistent both 
with the coupled Maxwell and matter-field equations and with the adiabatic 
regularization scheme. This is a nonperturbative calculation, which enables 
us to investigate the dynamical evolution of the interacting system of matter 
and the electromagnetic field. Careful computational work was essential in 
order to achieve numerical stability for this system since it consists of three 
time scales which are related to frequencies of the order of the mass of the 
produced particles (short time scale), the plasma oscillations of the electric 
field and current (medium time scale), and the degradation of the electric 
field to small values, by which point particle production is negligible (long 
time scale). It requires fine grids in momentum space and small time steps 
to solve the differential equations. 

Physical features like plasma oscillations and the plateau in the first pe- 
riod of the current evolution are understood by solving a classical relativistic 
system of particles in the presence of an electric field. From the value of the 
plasma oscillation frequency or the height of the first plateau of the current 
we obtain two separate estimates of the number of particles at the corre- 
sponding times. 

There is agreement between QED and a simple phenomenological model 
based on classical Boltzmann-Vlasov equations supplemented with a Schwin- 
ger WKB formula. This agreement is significantly improved when the source 
term is multiplied by a Pauli-blocking or Bose-enhancement factor. The 
phenomenological model enables us to obtain immediate physical insight and 
allows us perform calculations that are much faster and easier than those of 
the full QED formulation. In addition in the phenomenological model we 
understand the meaning of the number of particles, whereas in QED this 
concept is ill-defined as long as the interaction takes place. The agreement 
between the number of particles in these two calculations then allows us to 
assess the quality of the definition of the number of particles in the field 
theory in the presence of interaction. 
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The existence of a phenomenological model, in this simplified situation, 
with physical and numerical content very close to that of the exact field- 
theory treatment opens the way for studies that include spatial dependence, 
real radiation, and, eventually, color degrees of freedom. At the same time, 
it suggests that further application of phenomenological transport equations 
for systems involving pair production is in order. 
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Appendix: Particle spectra 

During the process of particle production, the particle number is not con- 
served, and we are not in an out-state region. However, it is possible to in- 
troduce an interpolating particle number operator, using the time-dependent 
creation and annihilation operators of the first-order adiabatic vacuum. This 
interpolating number operator has the property that if at t — the initial 
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state is equal to the first-order adiabatic vacuum state, then the number op- 
erator starts at zero. At late times, the first-order adiabatic number operator 
approaches the usual out-state number operator. 

The wave functions of the first-order adiabatic expansion, 

^(x,t) = e*~ft°\t), 



k 

i r*w k (t')<it' -iv m (t) 

(Q) e Jo kl > e iy k (t) 

Jk \ t ) 



(2o; k (t))V2 (2c k (t))V2' 
k (o) (x,t) = ft\t)e- ik -, (A.l) 

form an alternative basis for expanding the quantum field <3>, so that we can 
write 

= / [dk] (a k (t)0 k o) (x,t) + 6j c (t)0 k (o) (x,t)) . (A.2) 

Previously we expressed the field in terms of time-independent creation and 
annihilation operators: 

$(t, z, x ± ) = J [dk] (a fc k (x, i) + 6 k k (x, *)) (A.3) 

and 

0k(x, t) = / k (t)e ik - x , 0*(x, t) = f\(t)e-^. (A.4) 
Because both / and obey the Wronskian condition, it follows that 

(0 k ( X ,t),0 k ,(x,t)) = (0 k ° ) (x,t),0g ) (x,t)) = (27r) 3 5 3 (k-k / ), 

(0 k (x,t),0 k ,(x,t)) = (0* ( V,^0k°W))=O, (A.5) 
where the inner product is defined as 

* x {""at ~ Wj ■ (A ' 6) 

Using this we find that 

a k = (0 k (x, t), 0(x, t)) , &£ = - fe(x, t).0(x, *).) (A.7) 



At large times t, the time variation in the time-dependent creation and 
annihilation operators connected with the first-order adiabatic wave functions 
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becomes small. In that regime (which is not quite the out regime) we can 
determine these operators through the relations 

a k (t) = (#(x,t),0(x,t)) , bl(t) = - (0* (O) (x,t),0(x,t)) . (A.8) 

From these we can explicitly evaluate the Bogolyubov transformation at large 
times. In general we have 

a k {t) = a(k, t)a k + /3*(k, *)&!, b k (t) = a(k, t)b k + /3*(k, f)a£. (A.9) 

At late times, these operators become time-independent, and a and (5 are 
given by 



(0), 



AM) = *U m f-^A 



-/k , 



(A.10) 



The time-dependent interpolating particle number operator is defined by 
N(k,t)(2nf5 3 (k-k') = (4(t)a k (t)>, (A.ll) 
where N + = N_ = N. Thus, 

JV(k, t) = N(k) \a(k, t)\ 2 + (1 + N{k)) |/3(k, t) | 2 2i?e{a(k, *)/3(k, t)F(k)}. 

For the case of iV = F = we obtain at late times 

1 



(A.12) 



N(k,t) = 



2 1 ( Q k UJ k 



4 \n k u k 



(A.13) 



We note that our initial conditions ensure that initially this interpolating 
number operator is zero. 
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